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ABSTRACT 

Properties of a transformation method which has been developed for 
solving fluid dynamic problems on general two-dimensional regions are 
discussed. These Include the error in the construction of the trans- 
formation and applications to mesh generation. An error and stability 
analysis for the numerlceil solution of a model parabolic problem is 
also presented. 
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I. Introduction 


In solving boundary value problems In tvo-^lnenslonal regions, the 
treatment of curved boundaries has hindered the implementation of finite 
difference methods. Some have overcome this obstacle by using finite 
element methods at the expense of greater progranaulng complexity and 
computation time. Another alternative la to transform the problem to 
a region vhere finite difference methods are more suitable. The expense 
of a desirable computational region Is often a more complicated equation, 
or system of equations, to be solved. 

This report will examine the transformations developed by Chu {.!] and 
Thompson et al. [5] and [6] for solving fluid flow problems. The aim 
here Is to examine the basic transformation method and its construction 
and not the application to any physical problem. In this way an Indivi- 
dual with a particular problem may be better equipped to ascertain if 
this method is well suited for the solution of his problem. The sta- 
bility analysis for a model parabolic problem Is presented in the same 
spirit. 

II. The Elliptic Systems 

Many of the desirable features of these mappings over others which 
can be constructed derive from the fact that the mapping is obtained 
as the solution of an elliptic system of partial differential equations. 
In particular, the mapping will be differentiable at Interior points 
of the region. 

Let D be a bounded region whose boundary, denoted by 8D, Is the union 
of a finite number of disjoint closed contours C^,...,C^. Now D is of 



connectivity n end e simply-connected region cen be constructed from 
D by removing n-1 curves connecting different boundary components. 

These curves vlll be called branch cuts. On removing the branch cuts» 
the resulting region can be mapped by a one-to-one continuous function 
onto the Interior of a rectangle. The mapping can be extended continuous- 
ly to the closure of D, D, such that the extension Is one-to-one except 
on the branch cuts vfaere the mapping Is tvo-to-one. A typical boundary 
correspondence for a multiply-connected region Is given In Figure 1. 

With the boundary correspondence specified, the following boundary 
value problem can be stated. Determine functions C D which sat- 

isfy the semlllnear system 


5 - f (e. n) 

n ■ g (5, n) on D 


( 1 ) 


and the boundary conditions 


? " 4> (x.y) 

n - (x,y) on 3D. (2) 

The functions C snd ri are solutions of <1) in the following sense. 
They are to be treated as branches of solutions of (1) defined on a 
Rlemann surface. Thus the endpoints of the branch cuts are given by 
the boundary correspondence, but the Interior points of the cuts are 
determined by the solution of (1). The questions of existence and 
uniqueness of solutions of (1) and (2) will not be dealt with here. 
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Semlllnear ayatems have been atudled extensively for regions In the 
plane. Since the regions of a Rlemann surface that appear In this re- 
port can be mapped conformally onto a region of the plane, many of the 

results on semlllnear systems are still valid. It will be assumed that 
3D and the functions f* g, and tj; possess sufficient smoothness so that 
a solution of (1) and (2) exists which Is differentiable at all but a 
finite number of boundary points. In addition, the functions f and g 
will be chosen so that the image of any (x,y) In D belongs to R . This 
Is true If f = g = 0 by the Maximum Principle and. In general, restric- 
tions on the sign of f and g outside of R are sufficient for D 

to map Into R. 

In order that a transformation method be applicable to the solution 
of systems of partial differential equations. It Is necessary that the 
Jacobian of the transformation be nonzero on D. This Is the case, as 
will be shown next, for harmonic mappings of simply and doubly-connected 
regions . 

Let D be the simply-connected region bounded by the closed con- 
tour C. Decompose C into four arcs K^, K2» ^3* ^4 having only 
endpoints in common. The ordering and orientation of the arcs Is In a 
counterclockwise fashion around C, Let a and b be real numbers with 
a < b. 

Theorem 1. If C Is a harmonic function on the simply-connected 
region D and ^ > a on C ”* b on K^» € Is Increasing on K2, and ^ 

is decreasing on K^, then the gradient VC 0 on D. 
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Proof I Tho proof Is boiod on the Argument Principle. Suppose 
VC ■ 0 at (x^, y^) - 2 ^ In D. Since D Is simply-connected, the har- 
monic function C will have a single-valued harmonic conjugate C*. 

The analytic function 

W(e) - C (x,y) + 1 C* (x.y) 
will have V* (s^) - 0. Let 

u(e) - W(z) - Viz^) . 

The function (i) has a zero of multiplicity at least two at z - z^. By 
the Argument Principle, the change in the argument of u around C 
(or a curve In D arbitrarily close to C In case oj has a zero on C) 
must be at least 4 it. This contradicts the boundary values of C since 
C(x,y) - y^) assumes the value zero only twice on C. 

If D Is mapped to a rectangular region R by harmonic functions 
C and n with K^, K 2 , K^, mapping to the edges of the rectangle 
In a one-to-one manner, then both VC and Vn will be nonzero on D 
by Theorem 1. The fact that the Jacobian Is nonvanishing for simply- 
connected regions was mentioned by Godunov and Frokopov [2]. The validity 
of their argument requires Theorem 1 although this Is not stated. There 
Is also the question as to whether their argument, which uses conformal 
mappings, can be generalized to multiply-connected regions. An alternate 
proof Is given here ^Ich can be used for simply and doubly-connected 
regions. It is clear that the same reasoning can be used to prove that 
transformations of many regions of higher connectivity have nonvanlshlng 
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jACoblami. However, no proof haa bean atteaptad for raglona of arbi- 
trary connectivity due to the nuaerous naye of aaalgnlng boundary 
correspondancea and branch cute* 

Theorem 2. If the alnply-connec Ced region D Is mapped to a 
rectangular region R by the harmonic functions K nnd ri. then the 
transfoimatlon has a nonvanishing Jacobian on D. 

Proof: Let z ■ (x . y ) be an arbitrary point of 0. Then VC 0 

o o o 

at z and the curve 
o 

{(x»y)U(x,y) « C(x^. y^)} 

detemlnes local orthogonal coordinates (s.n) where s is the arc length 

parameter and n is the normal In the direction of decreasing C« In 

terms of these local variables, the Jacobian at any point of D is 

-C n • Now C and n are harmonic on D and £ < 0 since V£ 0. 

n s n s n 

The function n assumes Its minimum and maximum values on two arcs of 

D, say K 2 and K^. On and K^. n Is an increaslrg function of 

8. Therefore, except at a finite number of points. ^ 0 on 3D which 

Implies n > 0 on D. This proves that the Jacobian Is positive in D. 

8 

Interchanging K 2 and would result In a negative value. 

Let D be a doubly-connected region interior to the contour 
and exterior to C 2 * Let n be the harmonic function on D which 

assumes the boundary values n ■ a on C 2 and n “ b on with a < b. 

Let the boundary values of £ be prescribed such that £ increases 

from some value c to a larger value d as and C 2 are circum- 

scribed in the counterclockwise direction beginning at points z.^ and 
Z 2 on Cj^ and C 2 , respectively. Let £ be harmonic except on a 
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branch cut fron to »2 thara la a Jmtp of d > c In tha valua 

of the function. 

Thaoren 3. The tranafomation of the doubly-connected region D 
onto the rectan^dlar region R has a nonvanlahlng Jacobian on D. 

Proof; The function r) la constant on and C 2 . A simple 

relation betuean n and tha conformal napping of 1> onto an annular 
region shovs that the gradient of n la nonzero In D (see Ohtsuka 
[3, pp. 44-47]). The rest of the proof la similar to the proof of 
Theorem 2. Local coordinates (stn) are introduced so that the Jacobian 
takes the form Now < 0 on D and Cg la a single-valued 

harmonic function on D which Is nonnegative on 3D. Thus ^ 0 
on D and “5,^^ ^ 

The tranafomatlcna of the simply and doubly-connected regions have 
nonvanlahlng Jacoblana which together with the prescribed boundary 
correspondence guarantees that the mapping of D to R Is one-to-one 
and onto except on the branch cuts. A direct proof that the harmonic 
mapping of the doubly-connected region is one-to-one and onto has been 
given by P. D. Lax (private communication). 

The Theorems In this Section are not true for all transformations 
generated by the system (1). Nevertheless, the scope of our investi- 
gation will be broadened to include solutions of (1) and (2) with the 
additional hypothesis that the Jacobian of the transformation Is non- 
zero on D. A first step In that direr.tlon Is a remark on the effect 
of the functions f and g on the transformation. If f Is increased 
from 0 to a positive function, C becomes subharraonlc . Its values 
In D are decreased resulting In a transformation with higher resolu- 
tion near those boundary components where C assumes Its maximum and a 
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movement of the branch cuts of 4 In the direction of Increasing 
If f Is decreased from 0 to a negative function, higher resolution 
occurs near boundary components where 4 assumes Its minimum and the 
branch cuts move In the direction of decreasing S imi lar statements 

hold for n and g. This concept has been refined by Thompson et al. 
[6j to produce transformations with high resolution In various subsets 
of the region. 

The Inverse transformation from the rectangular region R to the 
region D Is also the solution of a system of partial differential 
equations. In fact, the system (1) is equivalent to the quaslllnear 
elliptic system 


a - 23 + Y x^^ + [f<5,n) x^ + g(5,n) x^] - 0 


a - 23 y^^ + Y y^^ + [f (C.n) y^ + g(C,n) y^3 - 0 


where 


2 2 
a « x„ + y 
n n 


V - 2 ^ 2 

Y-Xj 


■ ■‘{''n ■ V£ 


(3) 


Boundary conditions are obtained from (2) In the form 
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X • *(C.n) 


y ■ *!'(?, n) for (x.y) on 3D. 

Note that no boundary conditions are Imposed on the Image of the branch 
cute. The values of the functions x and y are set equal on any two 
subsets of 3R corresponding to the same branch cut. there Is 

a Dlrlchlet condition on part of 3R and a periodicity condition on 
the remainder. 

The first application of solutions of (3) was in the construction 
of Irregular curvilinear meshes. The system (3), with f = g s 0, 
was solved numerically by Winslow [7] to create a mesh for finite 
difference calculations. If a square mesh is constructed on R, the 
Image in D will be a curvilinear mesh. All boundary components 
and branch cuts lie on mesh lines. The curves of the mesh, on tdiich 
either C or n Is constant, Is said to generate a curvilinear 
coordinate system for D. The following list contains various pro- 
perties of the curvilinear coordinate system and Indicates how they 
are related to the coefficients In (3). 

A. Curvature of ^ = constant: a ' |x^ v„ - x^y^^i 

' nn n n nn' 

- 

C. Angle of Intersection of ^ constant and n » constant: Arctan 

D. Metric along ^ = constant: I dn| 

E. Metric along r| constant: |d?l 

When creating a mesh for either finite difference or finite element 


B. Curvature of t] = constant: Y 
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calculations, the values In C, D, and E often have an effect on the 
accuracy of the calculations. 

III. Mr.^erlcal Transformation 

In the numerical construction of the transformation, the system (3) 

Is discretized and solved on a rectangular or triangular mesh using non- 
linear SOR. Assuming that one desired to solve some partial differen- 
tial equation, or system of equations, on D, those equations would 
also be written In terms of the variables x and y. When the equa- 
tions (1) are used In this conversion, which has been the case In 
previous applications, any error In tna iumerical solution of (3) 
could produce additional error In the numerical solution of the trans- 
formed equations. E;iamples will now be presented where the error In the 
construction of tne transformation can be analyzed. 

Consider the harmonic mapping of the unit disc 

D • {(x»y) 1 x^ + < 1 } 

onto the interior of the square 

R «= 1 1 < 5 < 21 , 1 < n < 21} . 

In order to construct the mapping, suppose the system (3) Is dis- 
cretized on a square mesh of width 1. The resulting difference 
equations are solved at the interior mesh points of R by SOR Iteration. 
The mesh points on 3R were assumed to map to equally spaced points 
on 3D. The generated curvilinear mesh Is exhibited In Figure 2. 
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At the mesh points of 0, the c: act solutions of the equations In (1), 
5 and n» would have Integer values 1, 21 whenever the solu- 
tion of (3) Is exact. For comparison, the values of and hence 

n by symmetry, were computed from the Poisson Integral Formula at 
the Interior mesh points obtained from the solution of the difference 
equations. The values of ^ beginning at the lower left point of 
the mesh are given In the Table. Due to symmetry, only values In 
the left half plane are listed. Note that the maximum absolute error 
occurs near the boundary points where the Jacobian vanishes; that Is, 
the points Tfhich map to the vertices of the square. The error given 
In Figure 2 is normalized by dl'^ldlng all function values by the 
width of the square to compare with the following example. 

For a second example, a conformal mapping Is constructed. J.et 
D be the annular region 

D ■ {(x,y) I 1 < x^ + y^ < e^^} . 

A conformal mapping of D onto 

R - {(C.n) I o<C<i, o<n<i} 

Is given by the equations 

?(x,y) = log (x^ + y^) 

1 y 

n(x,y) - ^ arctan {^) . 

The transformation was constructed numerically as In the previous 
example resulting In the mesh of Figure 3. The coordinates of each 
mesh point in D were substituted in (^) and these values were com- 
pared with the coordinates of the mesh points on R. It was found 
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that the greatest error was In the construction of the function ri 
which has a jump across the branch cut. In all cases the constructed 
values of ^ at the mesh points of D were accurate to five signi- 
ficant digits. 

It Is hoped that these examples will be of some value to anyone 
planning to use this transformation method to solve a certain problem. 

At least some estimate of the magnitude of error In the construction 
of the transformation can be conjectured as well as the effect of branch 
cuts and boundary points where the Jacobian vanishes. 


IV. Stability and Discretization Error 

The solution of a time dependent system of partial differential 
equations by an explicit finite difference scheme requires a sta- 
bility restriction on the size of the time step. If the system of 
equations is transformed to a rectangular region, the stability analy- 
sis must be carried out on the finite difference analogs of the 
transformed equations. This will be illustrated by outlining the 
von Neumann stability analysis of the linearized vorticlty trans- 
port equation with the forward time - central space differencing 
scheme. Except for the transformation aspect, the following remarks 
parallel those in Roache [4, pp. 51-53]. Boundary conditions are 
neglected. 

Let CO be a function 6t 'x, y, and t which is a solution 
of the partial differential equation 


w^ = -uco - vto +uVco 
t X y 


(5) 
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for t > 0 and (x,y) In D where U|V, and y are constants. Under 
the transformation (1) • the equation (5) becomes 

“"t • -jr ‘*‘55 " ^ ^ “"en + ^ 

for <C»n) on R where 

U"uy„-vx„ + yJf 
n U 

V"vx^-uy^ + yjg 

and a, 3» Y» and J are as defined In (3). Suppose a first order 

forward difference on a mesh of width At Is used for the time derlva- 

tlve and second order central differencing on a square mesh of width 

h Is used for the spatial derivatives. The value of (d at the 

point (jh, kh) at time step n Is denoted by h- ’ The valae 

J » K 

at time step n+1 would be given by the difference equation as 


n+1 n At r„ / n n \ 

“ “j.k " 2hJ ^““j+l.k " J~l,k^ ^ 


, n n . 1 

^“j,k+l " ‘^J.k-l^^ 


Vl,k ■ ^ ‘^j,k^ " ^Vl,k+1 **^j-l,k-l 


~ Vl.k-1 " ^ ^‘^j.k+l ‘^j.k-l ■ ^ ‘^j.k^^* 


The functions U, V, a, 3, y are evaluated at (jh, kh) with the 
derivatives replaced by the appropriate difference expressions. 
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Applying the standard von Keuaann Fourier analysis to the above 
difference eque.tlon, necessary stability criteria are 


.. . J h 

^ iupnvi • 

and 

,2 i_2 

- 2 jy| (a + Y) 

Note that unlike the usual stability criteria in Roache [4]» the value 
of V* In (^) appears In both stability conditions since It appears 
In the expressions for U and V. Also, the quantities on the 
right of the inequalities depend on the point (jh, kh) and the In- 
equalities should be satisfied at each mesh point If the same size 
time step Is \ieed throughout R* In the case D * R and the trans- 
formation Is the Identity mapping, the inequalities reduce to the 
usual conditions for stability of the forward time - central space 
difference method of solving (5). 

A few straightforward observations will now be made concerning 
the formal discretization error In solving (5). In the expressions 
for the a’v'lvatlves 

1 / X 

” T <5'n “e ■ H 
“y * -r <*E “n ■ ’‘i • 

2 

the difference approximation is of order h If there are positive 
constants K, H, and N such that 
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( 6 ) 


|J( >K 
a < 

Y < NJ^ 

for all h > 0. Since < aY, the Inequalities (6) also Imply that 
the difference approximation for 

A - -4 (a - 2 + g «>„ 

2 

Is of order h . Therefore, the above method for solving the linearized 
vortlclty transport equation (5) Is first order In time and second order 
In space whenever (6) holds* 
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Table. I^xact values of C for computed difference values c£ (€»n). 



234 567 89 10 11 j 

2 

1.9? l_,95 l.97_ 1.98 1.9.9 1#_99 lt99 1.951 1.99 1.99 

3 

?t9.9 ?.9<H ?.99 ?«97 ?.98 ?.98 2.99 2.99 2.99 2.99 

4 

3.97 3.95 3,95 3.9fi 3,97 3,98 3.98 3.98_ _3,99. 3.9^ 

1 5 

4.98 4.97 ii.98 4.9t 4*97. .4..9.8 4,98 4,.98_ 4.98 4^98 

6 

5.99_ 5,98 5,57.- 5,97 5,9.7 5,98.-5,98 5.98 5,98 5,98 

7 

6.99- 6,99 8.91 _ 6,-98. -6, 98 6*98- 6,98 .6,98 _6.99 6.J93 

8 

8.00 7.99 7.99 7.98 7.98 7.99 7.99 7.99 7.99 7.99 

9„ 

9,0 0 8,99 8,99 8.99 8,99 8 . 99 8. 99 8,99_ 1,99 8,99 

IP 

10.00 10.00 10.80 10.00 9.99 10.00 10,00 1Q,.I|0 10,00 10,00 

11 

11,18. AjU.80_U, 0-0- 11,80 11.10 11,08-11,00 11,00 11.80 -11.00 

12 

12,00 12.00 12,00 12,00 12,01 12, 08 12,00 12.80 12,00 .12,00 

13 

13.00 13.01 13.01 13.01 13,01 13,01 13,01 13,01 13.01 13.01 

. M. 

14. C_0 .14,. 01 14,01. 1.4j,8? -14dl2- t4,ftJL 14,01 14, 0 1. 14.8 1 14,8 1 

15 

15,01 15,81 15,82 15,02 15,02 15,02 15,02 15.02 15,01 15.01 

16 

16,01 16,02 16,03 16.03 16.03 16.02 16.02 16,02 16*02 16,02 

17 

17.02 17,03 17.04 17.03 17.03 17,02 17.02 17.02 17.02 17*02 

18 

18,03 18.85 18,05. 18. 04 18. Q3 16.02 18.02 18.02 16.01 18,01 

19 

19ii.05_-19-,ft6 19 jl05. 19,.03 19,8219-, 02-19. 01 19.01 19,01. 19,81 

2Q 

20.08 28.01 28.03 20,1? 20,11 2t,ll 28,81 20, 81. 28,8 1-18,81^ 
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Figure 1* Transformation of a multiply-connected region onto a rectangular 
region. 
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Mesh size Error 

21 X 21 0.004 



Figure 2. Curvilinear mesh on a disc and maximum absolute error In 
mapping to a square of unit width. 
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Mesh size Error 

21 z 21 0.01252 

42 X 42 0.00290 

84 X 84 0.00067 



Figure 3. Curvilinear mesh on an annular region and maximum absolute 
error In mapping to a square of unit width. 
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